load motor.mat;
XtestNorm = motorNormalize(Xtest);
XtrainNorm = motorNormalize(Xtrain);

alpha = logspace(-8,0,100);
alpha = alpha';
M = 100;
N = size(XtrainNorm,1);
marginalSigma = 20;
phiPoly = ones(N,M+1);
phiRadial = ones(N,M+1);
logLikelyHoodPoly = zeros(size(alpha,1),1);
logLikelyHoodRadial = zeros(size(alpha,1),1);


%polynomial precalulations
for rowDesign=1:N
    for colDesign=2:M+1
        phiPoly(rowDesign,colDesign)=XtrainNorm(rowDesign,1).^(colDesign-1);
    end
end

wPoly = phiPoly'*phiPoly \ phiPoly' * Ytrain;

summationPoly = 0;
for i=1:N
    summationPoly = summationPoly + (Ytrain(i) - wPoly' * phiPoly(i,:)')^2;
end

summationPoly = summationPoly / (2*marginalSigma^2);

%radial precalulations
mu = linspace(-1,1,M);
sigma = 3*(mu(2)-mu(1));
for x=1:N
    for j = 2:M+1
        phiRadial(x,j) = exp(-(XtrainNorm(x)-mu(j-1))^2 / 2*sigma^2);
    end
end

wRadial = phiRadial'*phiRadial \ phiRadial' * Ytrain;

summationRadial = 0;
for i=1:N
    summationRadial = summationRadial + (Ytrain(i) - wRadial' * phiRadial(i,:)')^2;
end

summationRadial = summationRadial / (2*marginalSigma^2);


for alphaIndex=1:size(alpha,1)
    bigSigmaPoly = inv(sigma^(-2)* (phiPoly') * phiPoly + alpha(alphaIndex)*eye(M+1));
    mPoly = sigma^(-2)*bigSigmaPoly*phiPoly'*Ytrain;
    
    logOfbigSigmaPoly = sum(log(eig(bigSigmaPoly)));
    
    
    logLikelyHoodPoly(alphaIndex)=(M+1)/2*log(alpha(alphaIndex)) - N/2 * log(2*pi*marginalSigma^2) +...
        0.5*logOfbigSigmaPoly -...
        0.5*alpha(alphaIndex)* (mPoly') *mPoly - ...
        summationPoly;

    
    
    
    bigSigmaRadial = inv(sigma^(-2)* (phiRadial') * phiRadial + alpha(alphaIndex)*eye(M+1));
    mRadial = sigma^(-2)*bigSigmaRadial*phiRadial'*Ytrain;
    
    logOfbigSigmaRadial = sum(log(eig(bigSigmaRadial)));
    
    
    logLikelyHoodRadial(alphaIndex)=(M+1)/2*log(alpha(alphaIndex)) - N/2 * log(2*pi*marginalSigma^2) +...
        0.5*logOfbigSigmaRadial -...
        0.5*alpha(alphaIndex)* (mRadial') *mRadial - ...
        summationRadial;
    
end
figure(281);
semilogx(alpha, logLikelyHoodPoly);
figure(282);
semilogx(alpha, logLikelyHoodRadial);
